## 

rm(list = ls())

## 

library(tidyverse)
library(rdrobust)
library(pbapply)

## Declare covars

covs <- c('turnout_party_2009_btw', 'soz_vers_beschaeftigte_share',
          'pop_density_km2', 'migration_out_share')

## Federal election data

bt <- readRDS('data/data_federal.rds') %>%
  filter(!is.na(treated)) 

## List of outcomes

outcomes <- c("turnout_party", 
              "agg_left_party", 
              "agg_center_party", 
              "agg_right_party")

## Source helper function to tidy rdrobust output

source("code/tidy_rd.R")

## Drop 2009

bt <- bt %>% 
  filter(year > 2012)

## First differences

diff_df <- pblapply(outcomes, function(o) {
  out <- bt %>%
    filter(year > 2012) %>%
    pivot_wider(values_from = o, names_from = 'year', id_cols = 'ags',
                names_prefix = 'o') %>%
    mutate(diff = o2017  - o2013) %>%  dplyr::select(ags, diff) 
  ## Rename
  colnames(out)[2] <- o
  
  ## Return this
  out
}) %>%
  reduce(left_join) %>%
  left_join(bt %>% dplyr::select(ags, pop_dec_09, applies_census,
                                 one_of(covs), state_id) %>%
              distinct(ags, .keep_all = T)) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>%
  filter(applies_census == 1)

## Donut hole width list

donut_list <- seq(25, 150, 25)

## Estimate for all outcomes / all donut hole widths

out_rd_opt <- pblapply(outcomes, function(o) {
  
  ## Iterate over donut hole width
  
  lapply(donut_list, function(d) {
    
    ## Iterate over whether B-W is excluded or not 
    
    lapply(c(T,F), function(exclude_state) {
      
      if (exclude_state) {
        subset_select <- !diff_df$state_id == '08' &
          !between(diff_df$runvar, -d, d)
      } else {
        subset_select <- rep(T, nrow(diff_df)) &
          !between(diff_df$runvar, -d, d)
      }
      
      ## Estimate

      out <- rdrobust(y = diff_df[subset_select, o] %>% pull(!!o), 
                      x = diff_df$runvar[subset_select], c = 0,
                      covs = diff_df[subset_select, covs])
      
      ## Tidy and return
      out2 <- out %>% 
        tidy_rd(se_nr = 3) %>% 
        mutate(outcome = o, 
               exclude_state = ifelse(exclude_state, 
                                      'State excluded', 'All states')) %>% 
        dplyr::rename(se = std.error, pval = p.value, 
                      bw = bw_left_h,
                      bw_bias = bw_left_b) %>% 
        mutate(conf.low90 = estimate - qnorm(0.95)*se,
               conf.high90 = estimate + qnorm(0.95)*se,
               donut_bw = d)
      
      ## 
      
      out2
    }) %>% reduce(rbind)
  }) %>% reduce(rbind)
  
}) %>%
  reduce(rbind)

## Prep for plotting
## Load

results_rdd <- out_rd_opt %>%
  mutate(outcome = dplyr::recode(outcome, 
                                 `agg_left_party` = 'Left-wing parties',
                                 `agg_center_party` = 'Centrist parties',
                                 `agg_right_party` = 'Right-wing parties',
                                 `turnout_party` = 'Turnout')) %>% 
  mutate(outcome = factor(outcome, levels = unique(outcome)[c(2:4, 1)])) %>% 
  mutate(exclude_state = recode(exclude_state, 
                                `All states` = 'Full sample',
                                `State excluded` = 'B-W excluded')) %>% 
  filter(!donut_bw>200)


# Figure A.10 : Donut hole RD ----

p1 <- ggplot(results_rdd %>%
               filter(str_detect(outcome, 'Left')), 
             aes(donut_bw, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted',
             color = 'grey40')  +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, size = 0.65) +
  geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90), 
                width = 0, size = 1.05) +
  geom_point(shape =21,
             fill = 'white', size = 2.5) +
  ylab('RD estimate (p.p.)') + theme_bw()+
  facet_grid(outcome ~ exclude_state) +
  scale_x_continuous(breaks = seq(25, 175, 25)) +
  xlab('Donut hole radius (inhabitants)')
p1

